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The study investigates detonations with multiple quasi-steady velocities that have been observed 
in the past in systems with multi-peaked thermicity, using Fickett’s detonation analogue. A steady 
state analysis of the travelling wave predicts multiple states, however, all but the one with the highest 
velocity develop a singularity after the sonic point. Simulations show singularities are associated 
with a shock wave which overtakes all sonic points, establishing a detonation travelling the highest 
of the predicted velocities. 

Under a certain parameter range, the steady-state detonation can have multiple sonic points and 
solutions. Embedded shocks can exist behind sonic points, where they link the weak and strong 
solutions. Sonic points whose characteristics do not diverge are found to be unstable, and to be 
the source of the embedded shocks. Numerical simulations show that these shocks are only quasi 
stable. This is believed to be due to the reaction rates having been chosen to be independent of 
hydrodynamics which permits shocks anywhere behind a sonic point. 


I. INTRODUCTION 

The decomposition of certain reactive materials can 
occur in two or more distinct steps, characterized by mul¬ 
tiple peaks in the thermicity (effective rate of energy re¬ 
lease). Nitromethane-air detonations [T] and other usual 
fuels using NOa, as oxidizer [5] give rise to such multi¬ 
ple reaction zone detonation structures. Thermo-nuclear 
fusion reactions also occur in sequential steps. Detona¬ 
tions in degenerate white dwarfs undergoing supernova 
explosions of the Type la have three sequential steps 
where carbon, oxygen, and silicon undergo fusion [3]. Hy¬ 
brid detonations are self-sustained detonation waves in a 
mixture of reactive gases with suspended reactive dust 
and display two sequential reaction zones in the detona¬ 
tion structure[4]-[6]. The gas phase reaction first proceeds 
without influence from the solid phase, other than en¬ 
ergy used to heat the particles, and momentum lost to 
the solid phase by entraining the particles with the gas 
flow. The solid phase reacts exothermically once it has 
absorbed sufficient energy. 

A common feature of multi-peaked thermicity[7j sys¬ 
tems is the presence of endothermic processes coupling 
the multiple reactions. These losses can be manifested by 
heat and momentum losses to confining tube walls, mass 
divergence, or curved geometries. Losses can also be in¬ 
trinsic to the system, as they are in hybrid detonations 
for example where particle heating and drag withdraw 
energy from the gas phase. Experiments and numeri¬ 
cal simulations in these hybrid systems have shown that 
the selection rules and detonation wave pressure profiles 
depend intimately on the kinetics of the reactions in ei- 
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ther phase, the amount of energy release and hydrody¬ 
namic resistance of the particles to the gas phase motion; 
Zhang’s recent review [S] provides the state-of-the-art. 

In systems with a simultaneous exothermic and en¬ 
dothermic process, the competition between energy ad¬ 
dition and loss dictate the structure of the self-sustained 
wave[n]. Self-sustained, steady travelling detonations 
with losses have a surface of zero net thermicity within 
their reaction zone, where energy release is balanced by 
losses. This surface is sonic with respect to the deto¬ 
nation front, in order to avoid a singularity in the solu¬ 
tion, and information beyond it is unable to reach the 
front. The detonation’s velocity is therefore only influ¬ 
enced by the accumulated energy release between the lead 
shock and this sonic surface. The conditions govern¬ 
ing the propagation velocity of detonations with losses 
are known as the “generalized Chapman-Jouget (CJ)” 
conditions and are treated at length in most detonation 
textbooks [7l fTO] . 

In systems with multiple thermicity peaks, it is plau¬ 
sible that the generalized CJ condition be met multiple 
times, and it is unclear what governs the speed of the lead 
shock. Previous investigations have examined the solu¬ 
tions of detonations with multiple exothermic reactions 
and simultaneous losses. 

Veyssiere and Khasainov HU [ID numerically studied 
steady hybrid detonations using two-phase reactive Euler 
equations. They found three steady propagation regimes. 
One regime was driven solely by the gas reaction while 
particles remained inert in the driving region. The sec¬ 
ond regime had both gas and particle reactions driving 
the detonation front, consequently travelling faster than 
previous regime. Finally, the third regime was steady 
only under certain parameters. At these specific param¬ 
eters, the detonation propagated at the same velocity as 
the first regime, however, they found a shock wave em¬ 
bedded within the reaction zone. When off-parameters 
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were used, the embedded shock eventually overtook the 
detonation front and the velocity increased to that of the 
second regime, where reactions from the gas and solid 
phases drove the detonation together. 

Bdzil et al. [13] investigated a two-reaction case with 
losses due to shock curvature, using asymptotic expan¬ 
sions of the reactive Euler equations. They found a range 
of curvatures (losses) where two quasi-steady detona¬ 
tion velocities were possible. Time-evolution of a blast- 
initiated detonation showed that, with certain initial con¬ 
ditions, the blast would decay such that the flow would 
adopt the slower detonation velocity, then in some cases 
abruptly transition to the fast solution. 

A multiplicity of steady velocity solutions for detona¬ 
tions with multi-peaked thermicity were found in these 
studies, some unstable. Internal shocks were found to 
be transient in some cases and steady in others. The 
selection rules for the detonation structure and velocity, 
however, remain unclear. This is partially because the 
origin and transient of the internal shocks is uncertain. 

This investigation aims to clarify the steady reaction 
zone structure and velocity of detonations with multi- 
peaked thermicity. Solutions and selection rules for sit¬ 
uations where two sonic planes propagate at different 
speeds will be explored. The case where the two sonic 
planes propagate at the same velocity is complex and 
poorly understood and will also be studied. Experiments 
and simulations in hybrid systems HiaiSllIldl] have 
suggested that a double reaction zone structure is possi¬ 
ble, with an embedded shock located between two sonic 
planes. It is presently unclear if such a shock is stable, 
if it depends on external losses, as modelled by previ¬ 
ous authors, or is dependent on other three-dimensional 
effects in the experiments. The wave structure of such 
double-structure detonations, their stability, and origin 
of embedded shocks will be studied analytically and nu¬ 
merically. 

This investigation begins with the detonation model 
introduced by Eickett in the early 1980s [HKISI, which 
takes the form of the reactive Burgers equation with a 
reaction source term. The model neglects the rear facing 
pressure waves of gas dynamics, hence significantly sim¬ 
plifying the mathematical complexity of the description. 
The model retains the important physics of reactive com¬ 
pressible flows and its complex dynamics, namely that 
pressure waves receive amplification, modulated by the 
local rate of energy release, and can form shocks m- 
Eickett has already demonstrated how the model can re¬ 
produce the complex steady structure of eigenvalue det¬ 
onations in the presence of one exothermic and one en¬ 
dothermic reaction m- A similar simplified model has 
also been used by Earia and Kasimov [18] to investigate 
the effect of losses on detonations and their stability. The 
present study adds a second exothermic reaction, *.e., a 
system with two peaks of thermicity and an endother¬ 
mic reaction. Using this simplified mathematical model, 
the structure of the steady state solution can be studied 
analytically. 


Fickett’s detonation model is presented in section jn] 
and is then solved analytically for the steady state in sec¬ 
tion [ml Integral curves of the steady solution are stud¬ 
ied, and a Hugoniot-Rayleigh analysis is performed in 
section IV Unsteady numerical simulations of the model 
are then presented; the method is described in section jV] 
and the results are shown in lYll Solutions, stability, and 
selection rules are then discussed in section [Vll| in light 
of the analytical and numerical results. The conclusions 
are summarized in section rVIIIl 


II. MODEL: REACTIVE BURGER’S EQUATION 


The model is based on the reactive form of Burger’s 
equation introduced earlier by Eickett [T51 [TB] . The 
model is an analogue of detonations; it neglects rear¬ 
facing pressure waves, simplifying the mathematical com¬ 
plexity of the system, while retaining the important 
physics of reactive compressible flows and their dynam¬ 
ics. The hydrodynamic model is 


dtp + d^p = 0 (I) 

where a: is a Lagrangian coordinate, t is time, p represents 
density, and p represents pressure. In this study, the 
equation of state 





( 2 ) 


is used, where are the reaction progress variables which 
range from zero (unreacted) to one (reacted). The con¬ 
stants Qi represent the heat release when positive, and 
losses when negative. Two sequential exothermic reac¬ 
tions (subscripts I and 2) and one loss (subscript 3) 
are considered with state-independent rates. The first 
exothermic reaction and the heat loss begin once shocked, 
while the second reaction begins upon completion of the 
first reaction. Simple depletion reaction rates are chosen 
for the exothermic reactions while a constant rate is used 
for the loss 


n =dtXi = fci(I - Xi)''\ 

r 2 =dtX 2 = k 2 {l - A 2 )'"U and (3) 

^3 =dtX3 = fca, 

where ki are scaling constants, and Vi represent the re¬ 
action orders. 

Initial conditions ahead of the wave are uniform, taking 
the values po = Po = Xip = 0 for simplicity and without 
loss of generality. 

The parameters used in this study are listed in table 
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TABLE I: Parameter set used in study 


ordinary differential equations 


fci Qi 

k2 

Q 2 

Q 2 V2 

ks Qs 



(1 CJ point) 

(2 CJ points) 


1 0.5 0.5 

0.2 

0.5 

0.478239821826059 0.5 

0.1 - 0.8 



FIG. 1: The structure of self-sustained detonations 


III. STEADY ANALYTICAL SOLUTION 


A travelling wave solution can be deduced when the 
hydrodynamic equation[^is written in characteristic form 


dp da; 

-=».long- = p 


( 4 ) 


where t is time, and thermicity is defined as 




( 5 ) 


Figure [2 illustrates the structure of the self-supported 
detonation wave The detonation structure consists 
of pressure waves originating from the back (left), trav¬ 
elling along characteristics da;/dt = p and amplifying 
according to the characteristic equation Q. The am¬ 
plification is given by the local evolution of the reacting 
field, described by the reaction rate equations along 
the particle paths x = constant. These pressure waves 
coalesce and sustain a steady moving lead shock with 
velocity D. 

In order to seek the steady structure of the travelling 
wave solution illustrated in figure[2 the spatial variable x 
is changed to a shock fixed frame C, = x — Dt (with origin 
at the shock) for a detonation travelling with constant 
velocity D. The hydrodynamic equation 0 and reaction 
rate equations (|^ can then be transformed into a set of 


dp 

Tc 

d\i 

He 

d\2 

dC 

dC 


1 (T 

D p-D' 

n 


r2 


and 


D' 


( 6 ) 


Integrating these equations yields the analytical results 
for the wave structure 



Ai —1 - - I)(Ci - C) + , (8) 

A2=1- (^(i^2-l)(C2-C) + iy’" , (9) 

A3=^(C3-C)) (10) 

with Ci indicating where the reactions begin, thus = 0, 
C 2 = Cl + C 3 = 0 in accordance with the 

model. 

The travelling wave solution is isolated from the back 
when the limiting characteristic travels at the same veloc¬ 
ity as the steady lead shock, i.e., when dx/dt = p = D; 
this is the sonic criterion. For this limiting character¬ 
istic to travel at constant speed, it also requires van¬ 
ishing thermicity from the density differential equation 
([^. Thus the generalized Chapman-Jouget condition, 
denoted with the subscript CJ, is fulfilled when 

p = D and cr = i = 0 (II) 

i 

The solution exhibits sonic points, with respect to the 
front, at anywhere p = D. The relation between the 
detonation speed and the reaction, evaluated along the 
limiting characteristic, is obtained from the analytical 
relation for density 0 when the sonic portion of the 
generalized CJ condition 0 is met, and is 

(12) 

i 

This expression illustrates that the detonation velocity 
for the solution with sonic points is given by the net 
energy evolved from the lead shock to the sonic plane. 

The second portion of the generalized CJ condition 
0, the balance of reaction rates, permits the reaction 
progress values along the limiting characteristic to be es¬ 
tablished. Since the first and second exothermic reactions 
are sequential, the balance of rates must occur between 
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either the first and loss (third) reaction, or the second 
reaction and the losses, i.e. ri = or r 2 = r^. Denoting 
these sonic points as A and B respectively, two solutions 
are obtained 


Aia — 1 — 



, A 2 A = 0, and 


Asa — ^(Cs - Cl) + 


ki{vi - 1 ) 



for the sonic point closest to the shock, and 



Aib =1, A,b = 1 - (-|^) ’ . 

(») 

fca / / 

^2(1^2 ~ 1) yV ^2Q2/ 


for the second; recall Ci = Ca = 0- The sonic point 
positions are found by substituting these results into the 
wave structure equations andfTol), 


Ca 


Cb 


Ci- 


Ci + 


D _ 

ki{vi-l) A:2(lA2-l)yV ^2(32/ 



The solution is now complete for the detonation speed 
and reaction zone structure in closed algebraic forms. 


IV. STEADY ANALYTICAL RESULTS 

To illustrate the type of solution obtained, consider a 
numerical example with parameters such that I?a < Db . 
Figurej^shows four families of integral curves for the pa¬ 
rameters listed in table[T]with one CJ point. For these pa¬ 
rameters, Da = 0.59397 and Db = 0.6, while the equilib¬ 
rium detonation speed is Deq = = 0.447. 

The integral curves begin at the shock {C, = 0) with 
a value of p given by the inert shock jump conditions 
in Burgers equation, p = 2D, and proceed towards the 
burned side. 

A steady shock velocity D is chosen, paying attention 
to its relation to Dq and Da. When D > Db, the in¬ 
tegral curve (top most curve in figure [2a| ) does not in¬ 
tersect any sonic point; this is the over-driven solution 
which requires that the rear boundary be maintained 
at the corresponding value. The evolution of p is non- 
monotonous. Initially (zone 1), the heat release is greater 
than the loss, i.e. riQi > r^Q^, and the net positive 
thermicity leads to a positive density (and pressure) gra¬ 
dient, owing to the amplification of forward-facing pres¬ 
sure waves. The first zero in density gradient corresponds 


to when riQi = (point A), the first zero in thermic¬ 
ity. Towards the left (zone 2), riQi < r^Qs, the density 
increases as pressure waves are attenuated. Once the 
second reaction begins and overcomes the endothermic 
processes (zone 3), 7 - 2(32 > 7 - 3 ( 33 , density once again de¬ 
creases towards the left. After the second point of vanish¬ 
ing thermicity (point B) obtained when 7 - 2(32 = fsQsy the 
losses overcome the second exothermic reaction (zone 4). 
The last segment (zone 5) corresponds to when the losses 
terminate and the second reaction eventually comes to 
equilibrium. 

For a detonation speed corresponding to the largest 
eigenvalue, he., D = Db in this case, a sonic point ap¬ 
pears at point B, through the balance of the second re¬ 
action and the losses. Note that this sonic point is a 
saddle point, and both weak and strong solutions can be 
attained in the back, depending on the rear boundary 
conditions. This is the typical behaviour of pathological 
detonations, and well discussed by Fickett [16]. 

When D = Da, a sonic point occurs at point A, where 
the first reaction balances the losses. Note that the inte¬ 
gral curves corresponding to this solution terminate at a 
singularity. In the context of an unsteady solution, this 
signifies that a shock wave will form at the rear, which 
will eventually catch up to the lead shock. The unsteady 
solution presented in section [Vl| illustrates this transient. 
It can thus be asserted that the smaller eigenvalue is not 
stable, but can be established as an intermediate tran¬ 
sient. 


For detonation speeds lower than both eigenvalues, a 
steady solution does not exist, owing to the singularity 
established in the reaction zone, signifying the presence 
of a strong compression wave in the unsteady case. The 
equilibrium solution, where the detonation velocity is de¬ 
termined by the total heat release, is thus not possible. 


Inspection of the sonic point conditions (equations 13 


and 141 for determining the detonation speed eigenvalues 
(121 reveals the possibility of Da > Db. In this case, the 
largest eigenvalue. Da corresponds to the singularity free 
solution. 


The possibility of two sonic points in a steady solution 
is also present for a select parameter range such that 
Da = Db = Dab • This solution corresponds to a sin¬ 
gle integral curve passing through two saddle points, as 
shown in figure [2b| Parameters required for two simulta¬ 
neous sonic points need to satisfy a single constraint in 
the arbitrary choice of thermal-kinetic parameters, with 
a solution in closed form. Table HI lists the modifica¬ 
tion, in the second heat release for example, in order 
to achieve this condition. The solution obtained also ad¬ 
mits steady shock waves anywhere in the reaction zone 
structure. Since the shock speed in the Burgers’ and 
Fickett’s model is simply given by the average of left and 
right states, such shocks will travel exactly at the leading 
shock strength. However, the requirement that the shock 
be forward facing (impossibility of expansion shocks) re¬ 
quires these shocks to be behind a sonic point, as only 
there can a jump from the weak to the strong solution. 
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(b) Two CJ points, Da ~ Db — Dab 


FIG. 2: Family of steady integral curves for a detonation travelling from left to right with parameters found in table|^ 


A. Hugoniot and Rayleigh line analysis 

The multiplicity of solutions can also be represented 
in p-p phase space analysed by traditional Hugoniot-type 
arguments umis]. Constructing a solution requires con¬ 
necting possible unburned and burned loci (Hugoniots) 
with possible integral curves (Rayleigh lines) intersecting 
the loci of zero-thermicity {e.g. points A and B shown 
in figure [2a| ) . Such a representation offers further insight 
into the solution. 

For the simple Fickett model [16] considered here, the 
Hugoniots are simply given by the equation of state 
The unreacted shock Hugoniot (Ai = A 2 = A 3 = 0), the 
equilibrium Hugoniot (Ai = A 2 = A 3 = 1), the Hugoniot 
corresponding to the eigenvalue A (Ai = Aia, A 2 = A 2 A, 
A 3 = Asa) and the Hugoniot corresponding to the eigen¬ 
value B (Ai = Aib, A 2 = A 2 B, A 3 = Asb), are shown in 
the p-p phase space of figure for the single sonic point 
example studied in the previous section, whose integral 
curves are shown in figure |2a| Rewriting equation [^ as 
d{p — pD)/d( = 0, with boundary conditions at the lead¬ 
ing shock {p = 2D^, p = 2D, and labelled as N in figure 
[^, they form the Rayleigh lines 

p = pD. (15) 

which are shown in figure for the eigenvalue B and the 
over-driven solution of figure [^ 

The over-driven integral curve {D > subscript 

OD in figure starts at the leading shock, point Nqd, 
reaches point Aqd on the eigenvalue Hugoniot A, and 
then Bod on the eigenvalue Hugoniot B, then returns to 
Sod on the strong branch of the equilibrium Hugoniot. 
Points Aod and Bod are the local minima in the corre¬ 


sponding integral curve in figure as they correspond to 
the locus of zero thermicity, as explained above. 

The eigenvalue B integral curve {D = D^, subscript 
CJ,B) shown in figure is tangent to Hugoniot B. It 
starts at point Ncj.b, then intersects the Hugoniot A at 
point Acj 3 (first zero-thermicity point), then proceeds 
to the sonic point Bcj,b on the eigenvalue Hugoniot B, 
shown in detail at the top-left of figurej^ Since this point 
is a saddle point, as shown in the integral curve, figure 

the solution can then reach either the strong solution 
Scj.B or the weak solution Wcj,b- 

The eigenvalue A integral curve (not shown in figure 
1^ is tangent to Hugoniot A and slightly lower than that 
of eigenvalue B for this case. Since Hugoniot A is lower 
than B, the corresponding velocity is also lower, seen by 
the Rayleigh line equation Since this Rayleigh line 
never intersects Hugoniot B, a second zero thermicity 
point cannot be established in the system. The deto¬ 
nation speed selection rule and possible steady solution 
can thus also be made using the Hugoniot analysis. A 
regular solution requires intersection of the Rayleigh line 
with both zero-thermicity solutions. This can only be 
achieved for detonation velocities equal to, or larger than, 
the largest eigenvalue. 

In the case where two simultaneous steady sonic points 
are possible, the Hugoniot curves A and B overlap. 


V. UNSTEADY NUMERICAL METHOD 

Unsteady numerical simulations were used to study 
the stability and transients of the model. The do¬ 
main was discritized with a uniform grid spacing of 
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FIG. 3: Rayleigh lines and Hugoniot curves on a p-p diagram with the leading shock, equilibrium, and the two 
zero-thermicity (A and B) loci for the single steady sonic point properties; two sets of points are shown: an 
over-driven detonation {D > Db > Da, shown with subscript OD), and a detonation travelling at the CJ velocity 
corresponding to a sonic point at point B only (D = Db, shown with subscript CJ,B); N represents the unreacted 
shock condition, S and W correspond to the strong and weak solutions respectively 


Ax = 1/500. Time step size was determined using 
the Courant-Friedrichs-Lewy (CFL) condition such that 
CFL = 0.5 = Ax/{At X max(p)), where max(p) is the 
maximum value of p in the entire domain. A reaction 
threshold p > 0.001 was set to prevent unshocked gas 
from reacting. 

The Riemann problem was solved at every cell inter¬ 
face using a first-order Godunov method as described by 
Glarke et al. m- Time evolution was done using a first- 
order upwind method. 

An adaptive domain was used where cells were dynam¬ 
ically added ahead of the shock, while cells in the rear 
were removed once the flow had equilibrated. 


VI. UNSTEADY NUMERICAL RESULTS 
A. Single CJ point 

An example of the unsteady evolution of piston- 
initiated detonations, for conditions where the second 
sonic surface is slightly faster than the first (be. Db > 
Da), is shown with snapshots of the flow profile in figure 
1 ^ and with flow characteristics in figure 

A piston with velocity p = 0.35 was chosen in order 
to have an unsupported detonation but terminate at a 
value slightly larger than the weak solution. The initial 
condition can be seen in figure |4a| 


The first reaction initially drives a shock (figure |4b| ) . 
An over-expansion occurs once the first reaction weakens, 
followed by a recompression at the start of the second re¬ 
action, seen in figure The recompression forms into 
a shock (figure |4d[ ) which weakens and falls behind as 
the front continues to accelerate (figure [^, and a com¬ 
pression wave is formed at the piston as the loss reaction 
terminates just before the end of the second heat release. 
The piston is now behind the end of the reaction zone. 
The variable p increases at the beginning of the second 
reaction zone until it reaches the sonic condition, seen in 
figure 1^ The detonation front now travels at a speed 
Da corresponding to the sonic surface A. Meanwhile, a 
shock that travels faster than sonic surface A is created 
at sonic surface B (figure [4^, eventually penetrating into 
the first reaction zone. A truly steady detonation is even¬ 
tually established, travelling with a speed corresponding 
to the largest predicted eigenvalue Db, and with a new 
corresponding sonic surface B (figure [4i]). 

B. Two simultaneously steady CJ points 

A modification in the second heat release was made 
as discussed in section IIVI in order to allow for two si¬ 
multaneous steady sonic points. The possible integral 
curves for this example connecting the quiescent gas to 
rear conditions are shown in figure |2b| when D = Dab ■ 

The piston-initiated detonation of the double sonic 
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(j) Legend; left axis for p, right axis for Ai, horizontal axis for (^ = x — DAt 


FIG. 4: Wave profiles for piston-initiated detonation for the single sonic point case where < Dg (CFL=0.5) 


point solution, with the same piston velocity of 0.35 as 
the single sonic point case, is qualitatively similar to the 
single sonic point case previously discussed. The steady 
state emerged with only a single sonic point, as opposed 
to the two predicted by the integral curves. This discrep¬ 
ancy is reviewed in the discussion. 

The case with two simultaneous sonic points was then 


studied using the predicted steady-state solutions as ini¬ 
tial conditions to the unsteady simulations. All steady 
states start (C = 0) in the strong regime since the deto¬ 
nation is lead by a shock. Two options are available to 
the flow at each sonic point, following either the strong 
or weak solution, making for four possible steady states. 
This number is increase by the possibility of shocks join- 


































(a) At early times 



FIG. 5: Characteristic digram for a piston-initiated detonation for the single sonic point case (CFL=0.8) 


ing the weak and strong solutions. 

The shock-free initial steady-state curves will be re¬ 
ferred to with an abbreviated notation. This notation 
will break the flow profile into three sections: the rear 
solution, the solution between the sonic points, and the 
solution behind the leading shock. A strong solution will 
be abbreviated with the letter S and a weak solution with 
the letter W. For example, WWS would refer to a det¬ 
onation whose rear boundary follows the weak solution, 
has the weak solution between sonic points, and neces¬ 
sarily has the strong solution behind the leading shock. 
The WWS initial condition can be seen in figure [6a| 

Simulations initiated with shock waves behaved in the 
same manner as those which evolved shocks on their own 
and are therefore excluded. 


1. WWS initial condition 


The first simulation was initiated with the WWS inte¬ 
gral curve described in the previous paragraph and seen 
in figure The exact WWS solution was imposed 
on the domain at t = 0, then was numerically evolved 
through time. This solution can be achieved naturally 
with an under-driven piston, similarly to the single CJ 
point, piston-initiated detonation seen in figure [4^ Fol¬ 
lowing initiation, a disturbance travels rearwards (figure 
6 b) from the valley at the boundary between zones 2 
and 3. A shock forms at point B following the arrival 
of the disturbance (figure [6^. The shock travels faster 
than sonic point A, and eventually catches up to it and 
the front. During this transit, the detonation is slightly 
under-driven, while the final result is slightly over-driven 


and has a single sonic point at B, appearing qualitatively 
similar to case D = D-q in figure 

The unsteady simulation can also be visualized using 
the characteristic diagram in figure Recall that the 
model only admits the forward-facing family of charac¬ 
teristics (in the absolute frame of reference). This family 
of characteristics is shown in figure in the CJ deto¬ 
nation frame of reference. Characteristics outside the 
detonation’s domain of influence {i.e. to the right of the 
front) are omitted. Two quasi-steady runs of diverging 
characteristics are initially seen at points A and B. The 
internal shock, where characteristics of the same fam¬ 
ily interesect which results in a two-valued discontinuity, 
form at early times t < 50 and reach the detonation 
front between t = 400 and t = 450. The right-most char¬ 
acteristic, the detonation front, initially slopes slightly 
towards the left, showing it is slightly under-driven, but 
then slopes to the right when the internal shock reaches 
the front and it becomes over-driven. 


2. SWS initial condition 

Next consider the simulation initiated with a SWS inte¬ 
gral curve shown in figure |8a| Once again, a disturbance 
appears at the valley between zones 2 and 3, and trav¬ 
els backwards (figure [8b|) until it reaches the point B. A 
shock is created between the two sonic points (figure [8^. 
The shock links the supersonic flow (weak solution) be¬ 
hind the sonic point A to the new subsonic state (strong 
solution) ahead of the point B. The internal shock catches 
up to the detonation front (figures to in the same 
manner as described in the WWS case. Following this, a 
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(c) A shock forms and travels forward once the disturbance 
reaches point B (t = 44.0282) 


(d) Detonation is slightly over-driven following the internal 
shock’s arrival at the front {t = 450.477) 


FIG. 6: Wave profiles for WWS steady-state-initiated detonation for the double sonic point case 


new shock is created downstream of sonic point B to ad¬ 
just the flow to match the rear boundary condition. This 
new shock wave slowly distances itself from the front. 

Again, the detonation is slightly under-driven before 
the internal shock arrives at the front, after which it is 
slightly over-driven. 

The characteristic diagram is shown in figure]^ Rear 
of the domain has been omitted as it is largely uniform 
like that of figure Characteristics at point A initially 
diverge while those at point B converge. The internal 
shock is formed and the characteristics at point B then 
all become right-facing. The internal shock reaches point 
A, and a rear-travelling shock is then formed at point B 
and the characteristics at point B begin to diverge. 


Both sonic points initially have right-facing character¬ 
istics, seen in figure El They begin to diverge at sonic 
point B with the formation of the rear shock. 


4- WSS initial condition 

The fourth possibility, the WSS initial condition, is not 
shown due to its simplicity. Point A drifts into the strong 
solution, as seen in the SSS case, and the wave travels 
in a manner similar to the final state of figure |6d[ and is 
again slightly over-driven. 


VII. DISCUSSION 


3. SSS initial condition 


The simulation initiated with the SSS integral curve 
shown in figure |10a| where the strong solution is main¬ 
tained throughout. 

Once initiated, points A and B become over-driven 


(figure 10bI, seemingly independent of each other. The 


rise of point B is temporary and it returns to sonicity then 
drives the local flow to become supersonic (hgure jlO^ . A 
shock wave behind sonic point B is formed to adjust the 
flow to match the rear boundary conditions. The shock 
falls to the rear (figure lOdI as it did in the SWS case, 
while the detonation proceeds with a single sonic point. 
The detonation is slightly over-driven immediately upon 
initiation. 


The piston initiation of a detonation with parameters 
allowing only for a single steady C J point agree well with 
the steady integral curves. As the detonation front ac¬ 
celerates from the piston face, it begins with a velocity 
D < Da < Db- The detonation velocity temporarily 
stabilizes when D = Da as the creation of sonic point A 
stops signals from the rear boundary from affecting the 
front. A shock forms behind point A, as predicted by 
the singularity in the integral curve, as sonic point B is 
created behind the shock. During the transient of the 
internal shock, both sonic points A and B coexist. The 
internal shock eventually reaches the detonation front, 
accelerating it to D = Db and the detonation continues 
steadily with the single sonic point predicted in the in¬ 
tegral curves. Changing the parameters to allow for two 
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digram for the double sonic point detonation initiated with the WWS steady-state solution 

shown in figure 
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(c) An internal shock forms and travels forward once the 
disturbance reaches point B (t = 51.6496) 


(d) The internal shock still travels forward as flow at point B 
becomes slightly subsonic {t — 84.6758) 
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(e) Points A and B are temporarily over-driven following the (f) Point B becomes sonic and a backwards-travelling internal 
internal shock’s arrival at the front {t = 350.576) shock is formed {t = 450.395) 


FIG. 8: Wave profiles for SWS steady-state-initiated detonation for the double sonic point case 


steady CJ points {D = Da = Db = Dab) results in 
the same sequence of events, which does not lead to the 
predicted integral curve seen in figure 

Unsteady simulations were then initiated with the pre¬ 
dicted steady integral curves to understand the discrep¬ 
ancy, however, none of the steady profiles have proved to 
be stable as they all finished with a CJ point B and non- 
sonic point A. The formation of internal shocks, when 
created, always occurred at a sonic point, indicating that 
these sonic points may be unstable. 


Sonic points will be labelled with a method congruent 
to that used above to name steady-state solutions: sonic 
points lying between two weak solutio ns ( e.ff. point B in 
the WWS-initiated simulations, figure 6a) will be called 
a sonic point of type WW. Likewise, a sonic point be¬ 
tween two strong solutions is of type SS, a sonic point 
with a strong upstream solution and weak downstream 
solution is of type WS, and a weak upstream and strong 
downstream solution is of type SW. The top frames of 


figure [T^ serve as illustrations. 

All simulations initiated with the steady integral 
curves resulted in the creation of internal shocks, with 
the exception of the WSS initial condition. Point B was 
always the source these shocks, except when of type WS. 
As for point A, it quickly became over-driven when of 
type SS, or remained until shocked when of type WS. 
It appears sonic points of type WS are stable, unless 
shocked, while the others are unstable. 

The method of characteristics was used near the two 
sonic points to assess their stability. Consider the flow in 
the vicinity of a sonic point at steady state illustrated in 
figure |12| The top frames consist of the flow property p, 
which can be interpreted as the signal speed in Pickett’s 
analogy. The middle frame shows the corresponding in¬ 
stantaneous characteristics, or information paths, which 
have a slope df/d^ = l/{p—D). The bottom frame shows 
the resulting signal speed p after a period of time. The 
limiting characteristic (vertical line in the middle frame), 

























































440 

420 

400 

380 

t 

360 

340 

320 

300 

140 

120 

100 

80 

t 

60 

40 

20 

0 


12 



digram for a SWS steady-state-initiated detonation for the double sonic point case shown in 
figure]^ time between 150 and 300 is omitted in favour of clarity 






























13 



(c) Point B becomes supersonic {t = 42.3029) 


(d) A rearward travelling shock is created at point B 
{t = 100.673) 


FIG. 10: Wave 


profiles for SSS steady-state-initiated detonation for the double sonic point case 


travels along the sonic point at the same speed as the det¬ 
onation front. For the strong solution (when p > H), in¬ 
formation from the rear (left) is able to reach the limiting 
characteristic (and detonation front) as it travels faster 
than the front (it is hence subsonic in the with respect 
to the detonation. For the weak solution [p < D), in¬ 
formation cannot keep up with the front and is therefore 
carried towards the left in the detonation frame, away 
from the limiting characteristic (supersonic). 


For a sonic point of type WW (figure 12a) the char¬ 
acteristic diagram shows that any perturbation down¬ 
stream of the sonic point will easily be dissipated to the 
rear as the information speed (away from the front) in¬ 
creases in that direction. On upstream side, perturba¬ 
tions will accumulate and eventually disturb the sonic 
point. If a perturbation immediately upstream of the 
sonic point forces the signal speed p to locally become 
larger than D (he. adopt the strong solution just ahead 
of the sonic point), a shock wave must be formed to link 
the supersonic {p < D, weak) upstream solution to the 
new subsonic solution {p > D, strong) just ahead of the 
sonic point. Flow immediately upstream of the sonic 
point now satisfies the strong solution while maintaining 
the weak solution downstream of the sonic point. The 
sonic point of type WW has created an upstream shock 
and transitioned to a sonic point of type WS as seen in 
the bottom frame of figure |12a| 

Likewise, a sonic point of type SS (figure [l2b| ) will have 
upstream perturbations dissipated while those down¬ 
stream will accumulate. If the perturbation causes p < 


H, a downstream shock is created, joining the new weak 
solution immediately downstream of the sonic point to 
the strong solution further to the rear. A sonic point of 
type WS now exists, connecting the strong upstream so¬ 
lution to the new weak solution immediately downstream 
of the sonic point. The sonic point of type WW therefore 
creates a downstream shock and transitions to a sonic of 
type WS. 

Sonic points of type SW, shown in figure |12c| are also 
unstable. Perturbations accumulate on either side. An 
upstream perturbation acts in the same way as for a type 
WW sonic point, and a downstream perturbation in the 
same manner as for a type SS sonic point. Two shocks 
are created, one upstream and the other downstream of 
the sonic point, and the sonic point therefore adopts the 
WS form. 

Finally, the characteristics around a type WS sonic 
point (figure [12^ show that perturbations should be eas¬ 
ily dissipated away on either side. This explains the sta¬ 
bility of this type of sonic point seen in the simulations. 

From this analysis, it is expected that only the sonic 
points of type WS, where the characteristic paths diverge 
on either side, are stable. The rest should transform to 
create shock(s) and a stable sonic point of type WS. This 
leads to a further consequence: in order for two sonic 
points to coexist in a stable manner, they must be sep¬ 
arated by a shock wave. The location of this shock be¬ 
tween the two sonic points is arbitrary in the framework 
of this model. 

In these simulations, internal shocks (between the sonic 
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FIG. 11: Characteristic digram for the double sonic point detonation initiated with the SSS steady-state solution 

shown in figure [l0| 
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eventually changes into 



(a) Type WW sonic point 
results in an upstream shock 



eventually changes into 



(b) Type SS sonic point 
results in a downstream 
shock 



eventually changes into 



eventually changes into 



(c) Type SW sonic point 
results in shocks on either 
side 


(d) Type WS sonic point 
remains stable 


FIG. 12: Illustration of flow in vicinity of a sonic point: initial signal speed p (top), initial instantaneous 
characteristics (middle), and resulting signal speed p (bottom) 


points) eventually reached sonic point A and accelerated 
the detonation front. This instability is ascribed to the 
feature of the model, whereby any position of the inter¬ 
nal shock placed between the two sonic points is possi¬ 
ble. The detonation velocity was also found to be under¬ 
driven [D < Dab) prior to the internal shock’s arrival at 
point A, and over-driven after. The over-drive of point 
A, seen in figures |4^ and |10d[ is due to the det¬ 

onation over-drive. This is believed to be due to numer¬ 
ical diffusion. Increasing the resolution had a number 
of effects: reduced difference between expected (CJ) and 
measured detonation speeds, reduced difference between 
under- and over-driven detonation speeds, increased ac¬ 
curacy of shock-jump conditions across front, reduced in¬ 
ternal shock speed, and sharper features {e.g. point A 
of type SS). It is expected that a type SS sonic point 
A would cast a shock wave backwards between the two 
sonic points and become of type WS in the absence of 
this issue. 


VIII. CONCLUSION 

Detonations with two sequential heat releases and a 
concurrent loss were studied using Pickett’s detonation 
analogue. The study aimed to clarify the steady reaction 
zone structure and velocity when the two sonic planes 


propagated with different speeds (the general case). It 
was found that a steady detonation must travel at the 
maximum velocity dictated by the net heat release be¬ 
tween the front and each zero-thermicty point. Singular¬ 
ities were found in the flow structure at lesser velocities, 
and the lack of sonic point at larger velocities open the 
detonation front to the effect of perturbations. This sup¬ 
ports other studies [TTHI^ which found transitions from 
the lower detonation velocities to higher ones, but not 
vice-versa. 

Piston-initiated unsteady numerical simulations 
showed that the lower detonation velocities were quasi¬ 
steady. The formation of a forward-travelling shock and 
its subsequent arrival at the front caused the detonation 
to accelerate to its higher velocity. This phenomenon is 
believed to be the cause of the detonation front’s sudden 
acceleration in the previously mentioned by Bdzil et 
al.[T5]. 

The second goal of this study was to determine the sta¬ 
bility of the internal shock sometimes present in double- 
structured detonations. This was done by studying a 
special case where both sonic planes dictated equal deto¬ 
nation speeds. The presence of one sonic point in the flow 
allows the existence of an internal shock behind it, and 
the second sonic point should protect the internal shock 
from downstream perturbations. Internal shocks were 
found to be unsteady, likely due to the sequential reac- 
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tion rates being independent from the flow field, which 
allows internal shocks to exist anywhere. It is not clear 
how reaction dependence on hydrodynamics could an¬ 
chor internal shocks, further work is needed to establish 
in what capacity internal shocks can be stable, as sug¬ 
gested by Veyssiere and Khasainov fTTl IT^ . 

Finally, the origin of the internal shock was sought. 
Perturbations were found to accumulate at any sonic 
point with characteristics that not diverge on both sides, 
and lead to the creation of internal shocks. As a con¬ 
sequence, two sonic points may only co-exist in a stable 
fashion if they are separated by a shock. This clarifies the 
roots of double front detonations seen in previous work 
usiiiaiTi], and is also the mechanism by which deto¬ 
nations transition from slower, quasi-steady velocities to 


faster ones. 
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